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ABSTRACT 

The structure of constitutive equations for nonlinear multiaxial behavior of 
transversely isotropic fiber reinforced metal matrix composites subjected to proportional 
loading has been investigated. Results from an experimental program have been 
combined with numerical simulations of the composite behavior for complex stress states 
to reveal the full structure of the equations. It was found that the nonlinear response can 
be described by a quadratic flow-potential, based on the polynomial stress invariants, 
together with a hardening rule that is dominated by two different hardening mechanisms. 


INTRODUCTION 

Useful constitutive equations are given by the simplest expression that can model an 
observed behavior with a desired accuracy. Various aspects are stressed in different 
applications. In the initial conceptual design stage convenience is more important than 
precise prediction and simple expressions are used for scoping calculations. However, 
more accurate descriptions are required in the later stages of the design to optimize and 
verify the performance of the structures. The latter type will be addressed here. Metal 
matrix composites (MMC’s) reinforced with continuous ceramic fibers are strong and stiff, 
with only a weak nonlinearity in the response after matrix yield, when loaded in the fiber 
direction. In transverse tension and inplane shear the matrix carries a substantial part of 
the loading and the composite response is highly nonlinear when the matrix exhibits a 
nonlinear response. This strong anisotropic behavior imposes special requirements on the 
structure of the constitutive equations. 

Hill’s anisotropic yield surface has been used by fitting the free constants to a limited 
set of tests [Griffin et al, 1981] and [Bahei-El-Did et al, 1989]. The form of the yield 
surface is far too restrictive and the calculated responses are only valid for a very limited 
range of the stress space. Other constitutive equations are based on micro-mechanical 
models with simple assumptions of the local stress state in the composite; cf. 
[Aboudi,1982], [Min, 1981], and [Dovark and Bahei-El-Din, 1982]. Despite the simplifying 



[Aboudi,1982], [Min, 1981], and [Dovark and Bahei-El-Din, 1982], Despite the simplifying 
assumptions the resulting constitutive equations are quite lengthy and are unable to 
describe all essential features of the mechanical behavior with a sufficient accuracy. 
Detailed micro-mechanical models that include all the features of the composite 
behavior, cf. [Jansson, 1991], require numerical models that are of the same complexity 
as the structural problem. This implies that they are unsuitable for the analysis of 
structural problems because of the massive computational power required. 

A continuum description of the composite behavior is a more powerful approach for 
the analysis of structural problems. Spencer [1988] formulated quite general constitutive 
equations for the case when the plastic deformation is inextensible in the fiber direction. 
Hashin et al [1974] assumed that the response in the fiber direction is elastic and that the 
transverse response is non-linear. Lou and Shapery [1971] assumed that the longitudinal 
strain in the matrix is governed by the fiber deformation, wheareas Robinson and Duffy 
[1990] assumed that the plastic deformation of the composite conserves volume. The in- 
plane response of a MMC has also been fitted to a quadratic flow-potential and an 
associated hardening law by Kenagat et al [1987]. Sun and Chen [1989] later noticed that 
this fit had too high of a longitudinal yield stress and effectively gave a linear elastic 
response in the longitudinal direction. 

Many of the previous continuum formulations are extensions of constitutive equations 
developed for metals and assumptions appropriate for metals have been included in the 
formulation of the constitutive equations for the composites. These assumptions are far 
too restrictive for MMC’s. In this paper the response of a detailed numerical model of 
a MMC will be used to provide guidance for the formulation of the constitutive 
equations. The model complements the experiments by providing responses for multiaxial 
stress states that would be difficult to achieve experimentally. It also has the advantage 
that the stress states can be selected so that the structure of the constitutive equations 
can be identified in a rational way. 


OBSERVED MECHANICAL BEHAVIOR 

The material used in this study is Du Pont’s FP/AL [Campion et al., 1978]. It consists 
of continuous alumina fiber with a fiber diameter of approximately 20 Mn in uniaxial lay- 
up. The fiber volume fraction is 55% and the composite features a ductile Al-Li matrix 
with a strong fiber matrix bond. Mechanical properties have been reported by Jansson 

[1991]. 

The longitudinal tensile curve, Fig. la, has an initial linear response up to a stress of 
approximately 250 MPa, thereafter the matrix yields after which the slope of the curve 
shows a small decrease. The contraction in the transverse direction, Fig. lb, is also linear 
up to matrix yield, when a small increase is observed. 


2 


The transverse tensile curve, Fig. 2a, has a linear response up to a stress of 75 MPa. 
The subsequent nonlinearity is much more pronounced than for the longitudinal loading. 
The ductility of 1% is low compared to the matrix ductility of 30%. The low ductility was 
attributed [Jansson, 1991] to high hydrostatic stresses that build up in the matrix during 
the plastic deformation and reduces the ductility of the matrix. The contraction in the 
other transverse direction. Fig. 2b, increases after matrix yield and tends to a slope of - 1. 
The contraction in the fiber direction, Fig. 2c, decreases after matrix yield and tends to 
a very low rate. 

Shear response for different shear loading are shown in Fig. 3. Inplane shear is given 
by the curves for a = 0° in Fig. 3, transverse shear is given by a - 90° , and loading with 
equal components of inplane and transverse shear are given by a = 45°. All the responses 
for the different loadings are very close and can be assumed to be identical. The 
difference is less than the variation observed between different samples. The failure strain 
in shear is close to the matrix ductility because high hydrostatic stresses do not build up 
[Jansson, 1991], 

MODEL OF COMPOSITE MATERIAL 

The fibers are randomly distributed in the transverse plane of the composite. In the 
analysis the fibers are assumed to be long parallel cylinders arranged in a hexagonal 
array, Fig 4a, with the fibers orientated in the 3-direction. A perfect bond is assumed 
between fiber and matrix. The hexagonal array is the periodical array with mechanical 
properties that have the closest symmetries to a transversely isotropic composite. They 
have identical transformation properties when the constituents are linear elastic. The 
hexagonal array has a deviation from transverse isotropy when the constituents have 
nonlinear responses, [Jansson, 1992]. The deviation is strongest for inplane shear. For the 
present volume fraction the shear loading t ]3 cannot be expected to resemble the 
behavior of a composite with randomly distributed fibers because it activates slip on 
planes that are unconstrained by the fibers. 

The theory of homogenization, cf. [Keller, 1976] and [Jansson, 1992], states that the 
macroscopic response of a periodical composite is governed by effective properties given 
by a boundary value problem on the unit cell when the wavelength of the macroscopic 
loading is long compared to the wavelength of the microstructure. The displacement field 
can then be written to a leading order as 

u i (a - )= <e i ->Xj+« j p (jc)+u i ° (0 


where <e ;j> is the global average strain giving rise to a constant displacement field on 
the unit cell, u° is an arbitrary constant displacement, and u^ is a periodic 
displacement field that is equal on opposite sides of the unit cell. The macroscopic stress 
can be calculated from the local stress distribution in the unit cell as: 
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(2) 


<CT ij > = 



The effective stress-strain response for the composite is determined by giving the 
average strain as input. The periodical displacement field and the local stress distribution 
are solved for and the effective stress is determined by use of eq. (2). All the considered 
stress states are symmetric or anti-symmetric with respect to the y, and y 2 axis, shown 
in Fig. 4a. The periodical displacement field has an inversion symmetry at the point 
(yj=h/3/2, y 2 = b/2). This implies that only an eighth of the indicated unit cell in Fig. 
4a need be analyzed. The fibers are assumed to be linear-elastic and the matrix response 
is given by an elasto-plastic J 2 power law hardening rule. Because of material nonlinearity 
and complex geometry the finite element method has been used to solve the problem of 
the unit cell. The mesh is shown in Fig. 4b. The numerical procedure and boundary 
conditions are fully described in [Jansson, 1992], 

It was found by Jansson, [1991] that the calculated responses agree well with the 
experimental observations for a matrix with an initial yield stress of 75 Mpa and a 
hardening exponent of 5. From Figs. 1-3 it can be deduced that the model exhibits all the 
features of the observed mechanical behavior of the composite. 

CONSTITUTIVE EQUATIONS FOR COMPOSITE 

General Structure 

The total strain for the composite is assumed to be decomposed into a linear elastic 
and a nonlinear plastic part according to 

defdel+de* ( 3 ) 


The elastic part is given as 


e 




(4) 


where the compliance tensor C^, has the symmetry properties of a transversely isotropic 
material, that is given by five independent constants. As mentioned earlier the inplane 
and transverse shear moduli are close and can be assumed to be equal and this reduces 
the number of constants to four. The plastic strain will be assumed to be given by a flow 
potential 


4 



(5) 



dX 


with dA, as a plastic multiplier. The remaining part of the paper will address the form 
of the flow potential F and the structure of the hardening rules that is required to 
exhibit the observed and computed behavior of the composite. 


Formulation of Flow Potential 

The yield surface is frequently used as a flow potential in plasticity theories and 
provides advantageous features such as stability and variational formulations. The shape 
of the initial yield surface for the composite was calculated by Jansson [1991]. 
Intersections of the yield surface and a plane orientated in the ct 33 direction and at an 
angle <p with the a n - direction is shown in Fig. 5 where is the initial yield stress of 
the matrix. Experimental results confirmed the general dimensions of the surface but the 
data points were not sufficiently close to give the detailed shape. The calculated surface 
is convex and can be described closely by a quadratic polynomial in the stress 
components for all stress states except close to o n = -o 21 for which it has a sharp ridge. 
A ridge can be modeled by letting the surface be homogeneous in a non-integer exponent 
[Hill, 1979]. For a material exhibiting isotropic work hardening, the yield surface expands 
with a self-similar shape and surfaces for a constant amount of dissipated plastic work 
have self-similar shapes. Numerically calculated surfaces in the o u -a 22 plane are shown 
in Fig. 6 for different amounts of dissipated plastic work (W p ). The surface with zero 
dissipated energy is the initial yield surface. The sharp corners of the initial yield surface 
disappear after a small amount of plastic dissipation. The subsequent shape is smooth 
and the expansion is nonuniform with faster growth in the o n = a 22 direction than in 
the o n = -o 22 direction. 


The transition regime between the surface with sharp corners to the smooth surface 
is short and does not justify the complexity that would be required to include this 
transition. Disregarding the corner, the surface can then be approximated by a quadratic 
polynomial in the stress components. The polynomial invariants of the stress tensor tor 
a transversely isotropic material, cf. Green and Adkins [1960], are 


° ll +a 22> 


1 / \2 2 
r 33 ; ~(p n -o 22 ) +a 

4 


12” <J 13 +a 23’ 



when the symmetry axis is orientated in the 3-direction and the symmetry o- = a - of the 
stress tensor has been used. This set of invariants is not unique but irreaucible. The 
present form is the most convenient for the formulation of constitutive equations because 
the different invariants can be directly associated with different loadings. Including terms 
up to second order gives the following polynomial in the stress invariants; 
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( 6 ) 
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Letting 


F-1=0 


(7) 


for initial yielding; defines K L as the longitudinal yield stress, K T /2 as the yield stress 
for transverse biaxial tension; K rs as the yield stress in transverse shear; K ls as the 
yield stress in inplane shear; and K lT to be a constant defining the interaction of 
combined longitudinal and transverse loading. Convexity of the surface requires the 
negative sign in front of the interaction term together with the condition 


, KK r 


( 8 ) 


The plastic strain components are normal to the flow potential and are given by eq. (5) 
where dA is the plastic multiplier whose magnitude is governed by the hardening 
function. 


The suitability of using the surface F as the flow potential and the structure of hardening 
rules will be investigated using the stress-strain responses for different loading paths. The 
stress-strain relationships are generated by the numerical model that provide plastic strain 
increments and current stresses as output. The loadings paths have been selected such 
that ratios of the yield stresses KQ , that define the flow potential, can be determined 
conveniently in terms of strain increments and current stresses. 


For example, the normality rule (5) gives the strain increments resulting from a tensile 
load applied in the 1-direction: 


de 




dX 


de& 


ii 


ii 


K 2 t 2 4 


dX 
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Eliminating dA. and solving for K TS /K T gives 


^rs_ 1 


de^del , 

de p n -de> 22 


( 9 ) 


in terms of the plastic strain increments. Similar expressions have been derived for the 
other loading paths used in this study and are given in appendix A. 

The general trends for the ratios of the calculated yield stresses are given in Figs. 7 - 
10. In the following it is convenient to separate the loadings into an axi-symmetric part 
(with respect to the fiber orientation) that is given by the first three terms in the flow 
potential (6) and shear part that is given by the last two terms in (6). For axi-symmetric 
loading the ratios remain constant after a short initial transition period. The ratio given 
by the initial yield surface agrees very well for K LT /K, , Fig. 7, and is somewhat higher 
than the calculated steady state value for K LT /K T , Fig. 8. 

The computations indicate that the interaction between inplane (a 13 , a 23 ) and 
transverse shear (a 12 ) is quadratic. This is also supported by the experimental 
observations given in Fig. 3. The other computed stress strain curves are given in Jansson 
[1992], A quadratic interaction can be expected because local stress states in the 
composite are orthogonal and interaction through the matrix yield condition is quadratic. 

Transverse tension is a combination of axi-symmetric biaxial tension (o u =a 22 ) and 
transverse shear (a 12 ) in the present formulation, cf. eq. (6) . This type of loading forms 
additional requirements on the evolution of the yield stresses. Ratios of the yield stresses 
defining the axi-symmetric part of the flow potential are given in Fig. 8 for transverse 
loading. The ratios attain a constant value after a short transition period and this is the 
same trend as for the axi-symmetric loading. Curiously the a 22 loading has a long 
transition period. A representative example of ratios of yield stress for axi-symmetric 
loading to yield stress for shear loading is given in Fig. 9. It indicates that the yield 
stresses for the axi-symmetric part initially increases much faster than the yield stress for 
shear loading. The ratios tends to an asymptotic value which is much lower than the one 
given by the initial yield surface. 

It can be deduced that the normality requires that the ratios of the yield stresses 
defining the flow potential for axi-symmetric loading remain constant during the loading 
and are closely given by the initial yield surface. The ratio of the yield stresses for shear 
loading to yield stresses for axi-symmetric loading decrease gradually during the plastic 
deformation to saturation values that are significantly lower than the values given by the 
initial yield surface. 
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Hardening Rule 

A frequently used hardening rule in metal plasticity is isotropic work hardening with 
a self-similar expansion of the yield surface. Applying a similar idea here the incremental 
plastic work is 

dW*=o x fle J 0°) 


The normality rule (5) gives relations of the incremental plastic strain components 
d€jj in terms of the stresses a ^ and the plastic multiplier dA. Inverting the relations to 
get expressions for the stresses in terms of the plastic strain components and substituting 
them into (10) gives the following expression for the incremental plastic work 


dWH 


4- 


*L*r 

Klt 
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del ) 2 
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(de p n +de p 22 )de ] m 
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■4(de p 2 ) : 






I dk 


( 11 ) 


The requirement for positive dissipation is the same as the convexity condition (8), 
which is fulfilled by the initial yield surface. For isotropic hardening and proportional 
loading the plastic work may also be expressed as a function of a scalar equivalent plastic 
strain that can be identified from (11) as 


(de p f- 
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( 12 ) 




For a material exhibiting isotropic work hardening values of the function F with 
constants Kj given by the initial yield surface and stresses from different loadings 
will fall on one master curve when plotted as a function of the dissipated plastic work or 
the equivalent plastic strain. In Fig. 10 the value of F has been plotted versus the 
equivalent plastic strain for different loadings. It can be deduced that the results fall into 
three different groups of curves. One group consist of axi-symmetric loadings with respect 
to the fiber that also has the strongest hardening. The second group is given by the 
transverse and inplane shear loadings that has a weak hardening. The third group is given 
by the transverse tension that is a combination of axi-symmetric and shear loading and 
has an intermediate hardening. The results also indicate that pure hydrostatic loading 
causes the same magnitude of plastic strain as longitudinal tension. This implies that 
incompressibility should not be used as a constraint on the constitutive equations when 
the nonlinearity of the longitudinal response is included. This also leaves the freedom to 
fit the experimental observation better. 

These observations suggest that the hardening function should be separated into two parts 
as 


dK v _ dK h dK l T 

V® If® If® 
^LT 


=H A de 


p 

Ae 


(13a) 


dK. _ n 

-r"s (13b) 


where the superscript o stands for the initial yield surface and the effective strains in the 
hardening functions are defined respectively as 


(<e) 2 = 


*r 

2 

(de p n +de p 2 ) 2 +(de p 3 ) 2 +l 

/c r 

2 K l 

2 

^LT 


(de p n +de p 22 )de 33 ( l^ a ) 
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(14b) 


(de ^) 2 = (de \ 2 ) 2+ - (de p n -de p 22 ) 2 + (de p n ) 2 + (de P 2J ) 2 
4 


The first equation (13a) gives the hardening for axi-symmetric loading and referring to 
Fig. 10 is assumed to be linear in the effective strain. It is dominated by the constraint 
induced by the stiff fiber and the results in Fig. 1 1 indicate that H A is large. The second 
equation (13b) gives the hardening for shear loading and is highly nonlinear because the 
fibers constraints the shear loadings moderately. As stated earlier the inplane and 
transverse shear curves are close and can be assumed to be equal giving K xs = K IS = K s . 

It is obvious from the previous results that the flow rule, eq.(5), and the flow 
potential (6) together with the hardening rules, eq. (13), give excellent results for pure 
axi-symmetric or shear loading. It remains, to evaluate how well they describe the 
behavior under combinations of shear and axi-symmetric loadings. The value of the flow 
potential F has been plotted in Fig. 11a for different transverse loadings by using the 
hardening functions (13). The hardening functions have been determined from the axi- 
symmetric and shear responses given in Fig. 10. An exact fit requires the flow potential 
F to attain a constant value of unity throughout the deformation. It can be deduced from 
Fig. 11a that the hardening is not sufficiently strong to keep the value of F constant and 
equal to one, especially for stress states close to transverse biaxial tension. 


This discrepancy can be overcome by introducing cross hardening terms in the 
hardening functions (13) with 


dK T _ dK L 


^ =H A de Ac +H Asde 

*LT 


p 

Se 


(15a) 


^H(el)de +H SA de ? Ae (15b) 


where the cross hardenings have been assumed to be linear in the effective strains for 
simplicity. A cross hardening of H SA = 0.2H A in the shear hardening causes the value 
of F to remain constant, Fig. lib. This effect can also be established by a cross-hardening 
of H as =30 H a in the axi-symmetric hardening, Fig. 11c. 

The most preferable of the two cross hardenings is the one that gives the closest 
trend to that required by the normality shown in Fig. 9. The history of Kg/Kp with no 
cross hardening is given in Fig. 12a. A cross-term in the shear hardening causes the ratio 
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to decrease more slowly, as indicated in Fig. 12b, while a cross-term in the axi-symmetric 
hardening causes the ratio to decrease faster, Fig 12c. A comparison of the results in 
Figs. 12a-c with Fig. 9 indicate that the axi-symmetric hardening has the most appreciable 
affect when a single cross hardening term is used. 

DETERMINATION OF CONSTANTS 


A longitudinal tensile test provides the initial yield stress K L ° and the interaction 
constant K LT ° is given through the transverse contraction as shown by the relation 


IS® — IS® 

a lt~ a l 


N 


2de% 


de 


11 


(16) 


derived from the normality condition. The constant II A in the axi-symmetric hardening 
function can thereafter be determined from the longitudinal tests as 


”a= 



33 


(17) 


The inplane shear tests provides the initial yield stress and the nonlinear shear 
hardening slope H s (e Se P )* 


The initial yield stress for transverse biaxial tension can now be determined from the 
transverse tensile test by use of the initial yield stress in shear as 



(18) 


where a T ° is the initial yield stress in transverse tension. The values of the initial yield 
stresses are restricted by the requirement of positive dissipation (8) and the initial yield 
stress for transverse biaxial tension is also higher than the initial yield stress in shear. 

The cross-hardening term H SA can be estimated from the transverse stress strain curve 
by satisfying the yield condition close to saturation. This calculation gives the relationship 
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( 19 ) 


where 


//< 


SA 
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-J H(e)de 
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( 20 ) 


The other cross-hardening term H AS has a strong effect on the contraction in the 
unloaded transverse direction and can be included when it is appreciable. 

The constitutive equations include many constants and a direct systematic 
determination of the constants cannot be expected to give the best possible description 
of the mechanical behavior. The response have to be evaluated and compared with the 
experimental data for the different loadings. It is likely that the description can be 
improved by adjusting some of the constants. Experience has shown that the best fit is 
usually obtained by using interactive computer graphics to display the results and a 
human link in the loop to make the final adjustments. 

CONCLUSIONS 

It was found that a quadratic flow potential based on the polynomial stress invariants 
can be used to describe the mechanical behavior of fiber reinforced composite subjected 
to proportional loading. This formulation requires that the hardening function is 
separated into two parts: a strong linear hardening for axi-symmetric loading, that is 
dominated by the stiff fiber, and a weak nonlinear hardening, that is dominated by the 
matrix, for the inplane and transverse shear loadings. Calculations and experiments 
indicate that the composite can be assumed to have identical responses for the two shear 
loadings. 

Cross-terms in the hardening functions improves the accuracy for transverse tension 
by increasing the hardening and rate of constraint in the fiber. The cross-hardening can 
be neglected if an error of 20 % in stress can be accepted in transverse tensile or shear 
response. 

The plastic strain components for longitudinal tension and hydrostatic loading are of 
the same order. Hence, a consistent simplification of the constitutive equations is to 
assume inextensibility in the longitudinal direction when the plastic deformation is 
assumed to be incompressible. This assumption makes the constitutive equations 
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particularly easy and suitable for hand calculations. It implies that the composite can only 
deform in inplane and transverse shear. 

In a subsequent paper these constitutive equations will be implemented into a general 
finite element program and the performance of the equations will be evaluated for 
representative structural components. 
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APPENDIX A 
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Transverse biaxial tension and longitudinal plane strain 
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Transverse biaxial tension 
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Longitudinal tension and transverse plane strain 


a 33 , de p n =de p 22 = 0 


( 40 ) 


de p u =de p 22 = 


, a il +a 22 a 33 


K r 


'<a=o 


(A.9) 


gives 




LT 


'33 


Kj ^ 4o 


li 


(A. 10) 


16 


Transverse tension 
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Transverse loading close to biaxial tension 
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Transverse tension, a) tensile curve, b) contraction in other transverse 
direction, c) contraction in longitudinal direction. 
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Fig. 3 Shear loading. Inplane shear a = 0, transverse shear a = n /2 , and a = 

tt /4 equal components of inplane and transverse shear. 
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Fig. 4 




Model of composite, a) Hexagonal array, b) Finite element mesh of unit 
cell. 
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Fig. 5 Initial yield surface in the cr n -a 22 -a 33 space. 
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Ratio of yield stresses for axi-symmetric part of flow-potential for axi 
symmetric loading. 
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Ratio of yield stresses for axi-symmetric part of flow-potential for axi 
symmetric and combinations of axi-symmetric and shear loadings. 
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